In-class Ex10

1 Modelling Spatial Interaction: Huff and MCI Models

1.1 Overview

Market area analysis defines and studies where customers originate, using demographics, spending, competition, and observation to guide location, marketing, and expansion decisions. The Huff model estimates store choice probabilities and market share from attractiveness and distance, helping businesses delineate trade areas and forecast revenue among competing outlets and potential profitability.

1.2 Getting Started

Before we getting started, it is important for us to install the necessary R packages and launch them into RStudio environment.

pacman::p_load(sf, tidyverse, MCI, tmap, ggstatsplot)

1.3 Importing geospatial data

stores <- st_read(
  dsn = "/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/In-Class_Ex10/data/geospatial",
  layer = "STORES_TC")

Reading layer `STORES_TC' from data source 
  `/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/In-Class_Ex10/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 28 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 202702.9 ymin: 2662196 xmax: 232773.2 ymax: 2693557
Projected CRS: TWD97 / TM2 zone 121

town <- st_read(
  dsn = "/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/In-Class_Ex10/data/geospatial",
  layer = "TOWN_TC")

Reading layer `TOWN_TC' from data source 
  `/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/In-Class_Ex10/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 28 features and 7 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 194797.1 ymin: 2654887 xmax: 240450.2 ymax: 2703838
Projected CRS: TWD97 / TM2 zone 121

village <- st_read(
  dsn = "/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/In-Class_Ex10/data/geospatial",
  layer = "VILLAGE_TC")

Reading layer `VILLAGE_TC' from data source 
  `/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/In-Class_Ex10/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 622 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 194797.1 ymin: 2654887 xmax: 240450.2 ymax: 2703838
Projected CRS: TWD97 / TM2 zone 121

The output print above shown that the three data sets are in TWD97 / TM2 zone 121 (i.e epsg: 3826) projected coordinates systems.

Important

For any geospatial analysis, it is important to check and confirm that all the geospatial data are in the similar projected coodinates system.

tmap_mode("view")

ℹ tmap mode set to "view".

tm_shape(town) +
  tm_polygons() +
tm_shape(village) +
  tm_polygons() + 
tm_shape(stores) +
  tm_dots()

Registered S3 method overwritten by 'jsonify':
  method     from    
  print.json jsonlite

tmap_mode("plot")

ℹ tmap mode set to "plot".

1.4 Importing aspatial data

Next, read_csv() of readr package will be used to import the two aspatial data sets into R environment.

pop2019 <- read_csv(
  "/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/In-Class_Ex10/data/aspatial/pop2019.csv")

Rows: 625 Columns: 11
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (5): COUNTY, TOWN, V_ID, VILLAGE, INFO_TIME
dbl (6): COUNTY_ID, TOWN_ID, H_CNT, P_CNT, M_CNT, F_CNT

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

store_sales <- read_csv(
  "/Users/cktan/Desktop/SMU/01_Geospatial Analytics (ISSS626)/Hands-on_Ex/In-Class_Ex10/data/aspatial/Store_sales_data.csv")

Rows: 5499 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Store Code, POD code
dbl (3): Bills, Amount, Ave Bill

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

summary(store_sales)

  Store Code          POD code             Bills            Amount       
 Length:5499        Length:5499        Min.   :   1.0   Min.   :      0  
 Class :character   Class :character   1st Qu.:  37.0   1st Qu.:  26702  
 Mode  :character   Mode  :character   Median : 186.0   Median : 139643  
                                       Mean   : 331.9   Mean   : 236787  
                                       3rd Qu.: 457.0   3rd Qu.: 331952  
                                       Max.   :7764.0   Max.   :6076060  
    Ave Bill      
 Min.   :    0.0  
 1st Qu.:  661.3  
 Median :  730.5  
 Mean   :  825.1  
 3rd Qu.:  862.2  
 Max.   :24356.0  

1.5 Data Preparation

1.5.1 Aggregating store sales

Code chunk below is used to aggregate the store sales data at Store Code.

store_sales_aggre <- store_sales %>%
  group_by(`Store Code`) %>%
  summarise(
    "Total_Bill" = sum(Bills, na.rm = TRUE),
    "Total_Sales" = sum(Amount, na.rm = TRUE))

1.5.2 Relational join

In the code chunk below, left_join() of dplyr package is used to append data from store_sales data.frame onto stores sf data.frame. At the same time a new field called Store_size has been created.

stores <- stores %>%
  left_join(store_sales_aggre,
            by = join_by(
              Store_CD == `Store Code`)) %>%
  mutate(Store_size = 100)

Next, summary() is used to check if the are any missing values in the data.frame.

summary(stores)

  Store_Name         Store_Addr          Store_CD           Total_Bill   
 Length:28          Length:28          Length:28          Min.   : 1767  
 Class :character   Class :character   Class :character   1st Qu.: 5698  
 Mode  :character   Mode  :character   Mode  :character   Median : 8634  
                                                          Mean   : 8021  
                                                          3rd Qu.:10349  
                                                          Max.   :14468  
  Total_Sales               geometry    Store_size 
 Min.   :1475963   POINT        :28   Min.   :100  
 1st Qu.:4044821   epsg:3826    : 0   1st Qu.:100  
 Median :6116336   +proj=tmer...: 0   Median :100  
 Mean   :5542852                      Mean   :100  
 3rd Qu.:6959911                      3rd Qu.:100  
 Max.   :8234234                      Max.   :100  

1.5.3 Tidying pop2019 data.frame

glimpse(pop2019)

Rows: 625
Columns: 11
$ COUNTY_ID <dbl> 66000, 66000, 66000, 66000, 66000, 66000, 66000, 66000, 6600…
$ COUNTY    <chr> "臺中市", "臺中市", "臺中市", "臺中市", "臺中市", "臺中市", "臺中市", "臺中市", "臺中市…
$ TOWN_ID   <dbl> 66000110, 66000110, 66000110, 66000110, 66000110, 66000110, …
$ TOWN      <chr> "大甲區", "大甲區", "大甲區", "大甲區", "大甲區", "大甲區", "大甲區", "大甲區", "大甲區…
$ V_ID      <chr> "66000110-002", "66000110-009", "66000110-021", "66000110-00…
$ VILLAGE   <chr> "大甲里", "中山里", "太白里", "孔門里", "文曲里", "文武里", "日南里", "平安里", "江南里…
$ H_CNT     <dbl> 454, 1781, 406, 383, 355, 2071, 2039, 785, 461, 559, 769, 80…
$ P_CNT     <dbl> 1255, 5252, 1471, 1001, 1264, 6687, 6739, 2642, 1618, 1977, …
$ M_CNT     <dbl> 626, 2533, 758, 538, 640, 3244, 3406, 1310, 838, 983, 1273, …
$ F_CNT     <dbl> 629, 2719, 713, 463, 624, 3443, 3333, 1332, 780, 994, 1239, …
$ INFO_TIME <chr> "108Y12M", "108Y12M", "108Y12M", "108Y12M", "108Y12M", "108Y…

glimpse(village)

Rows: 622
Columns: 11
$ VILLCODE   <chr> "66000110029", "66000120001", "66000120002", "66000120003",…
$ COUNTYNAME <chr> "臺中市", "臺中市", "臺中市", "臺中市", "臺中市", "臺中市", "臺中市", "臺中市", "臺中…
$ TOWNNAME   <chr> "大甲區", "清水區", "清水區", "清水區", "清水區", "清水區", "清水區", "清水區", "清水…
$ VILLNAME   <chr> "建興里", "鰲峰里", "靈泉里", "清水里", "文昌里", "南寧里", "西寧里", "北寧里", "中興…
$ VILLENG    <chr> "Jianxing Vil.", "Aofeng Vil.", "Lingquan Vil.", "Qingshui …
$ COUNTYID   <chr> "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B",…
$ COUNTYCODE <chr> "66000", "66000", "66000", "66000", "66000", "66000", "6600…
$ TOWNID     <chr> "B11", "B12", "B12", "B12", "B12", "B12", "B12", "B12", "B1…
$ TOWNCODE   <chr> "66000110", "66000120", "66000120", "66000120", "66000120",…
$ NOTE       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ geometry   <MULTIPOLYGON [m]> MULTIPOLYGON (((211953.9 27..., MULTIPOLYGON (…

The plot above shown that the village ID (i.e. V_ID) in pop2019 data.frame is not in the same structure as the village code in village sf data.frame.

Code chunk below will be used to remove the “-” sign from the values.

pop2019 <- pop2019 %>%
  mutate(V_ID = str_remove_all(
    V_ID, "-"))

Now, we are ready to append data in pop2019 data.frame onto village sf data.frame by using left_join() of dplyr package.

village <- village %>%
  left_join(pop2019,
            by = join_by(
              VILLCODE == V_ID))

1.5.4 Check for missing or zero values

summary(village)

   VILLCODE          COUNTYNAME          TOWNNAME           VILLNAME        
 Length:622         Length:622         Length:622         Length:622        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
   VILLENG            COUNTYID          COUNTYCODE           TOWNID         
 Length:622         Length:622         Length:622         Length:622        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
                                                                            
   TOWNCODE             NOTE             COUNTY_ID        COUNTY         
 Length:622         Length:622         Min.   :66000   Length:622        
 Class :character   Class :character   1st Qu.:66000   Class :character  
 Mode  :character   Mode  :character   Median :66000   Mode  :character  
                                       Mean   :66000                     
                                       3rd Qu.:66000                     
                                       Max.   :66000                     
                                       NA's   :5                         
    TOWN_ID            TOWN             VILLAGE              H_CNT     
 Min.   :6.6e+07   Length:622         Length:622         Min.   : 180  
 1st Qu.:6.6e+07   Class :character   Class :character   1st Qu.: 724  
 Median :6.6e+07   Mode  :character   Mode  :character   Median :1335  
 Mean   :6.6e+07                                         Mean   :1592  
 3rd Qu.:6.6e+07                                         3rd Qu.:2195  
 Max.   :6.6e+07                                         Max.   :7599  
 NA's   :5                                               NA's   :5     
     P_CNT           M_CNT          F_CNT       INFO_TIME        
 Min.   :  487   Min.   : 257   Min.   : 230   Length:622        
 1st Qu.: 2139   1st Qu.:1095   1st Qu.:1053   Class :character  
 Median : 3862   Median :1911   Median :1934   Mode  :character  
 Mean   : 4545   Mean   :2234   Mean   :2311                     
 3rd Qu.: 6151   3rd Qu.:3022   3rd Qu.:3171                     
 Max.   :18652   Max.   :8741   Max.   :9911                     
 NA's   :5       NA's   :5      NA's   :5                        
          geometry  
 MULTIPOLYGON :622  
 epsg:3826    :  0  
 +proj=tmer...:  0  
                    
                    
                    
                    

The output print above revealed that the are five missing values in H_CNT, P_CNT, M_CNT and F_CNT fields.

Hence, code chunk below is used to remove the five records with missing values from village sf data.frame.

village <- village %>%
  filter(!is.na(P_CNT)) 

any(is.na(village$H_CNT))

[1] FALSE

any(is.na(village$P_CNT))

[1] FALSE

any(is.na(village$M_CNT))

[1] FALSE

any(is.na(village$F_CNT))

[1] FALSE

1.6 Computing Distance Matrix

We will calculate the centroid of each village using the code chunk below. Then, compute the distance between the centroids of the villages and the stores.

village_center <- st_point_on_surface(village)

Warning: st_point_on_surface assumes attributes are constant over geometries

distmat <- st_distance(village_center, stores) 

Next, the code chunk below is used to convert the distance matrix into long data.frame format

rownames(distmat) <- village$`VILLCODE`
colnames(distmat) <- stores$`Store_CD`

distmat_long <- distmat %>%
  as.data.frame() %>%
  rownames_to_column(var = "VILLCODE") %>%
  pivot_longer(
    cols = -VILLCODE,
    names_to = "Store_CD",
    values_to = "distance"
  ) %>%
  mutate(distance = as.numeric(distance)) 

1.7 Preparing Interaction Matrix

Code chunk below will be used to build the interaction matrix required by MCI package.

interaction_matrix <- distmat_long %>%
  left_join(
    stores %>% 
      st_drop_geometry() %>% 
      select(Store_CD, Store_size),
    by = "Store_CD") %>%
  relocate(
    distance, .after = last_col())

summary(interaction_matrix)

   VILLCODE           Store_CD           Store_size     distance      
 Length:17276       Length:17276       Min.   :100   Min.   :  108.5  
 Class :character   Class :character   1st Qu.:100   1st Qu.: 6678.7  
 Mode  :character   Mode  :character   Median :100   Median :11965.3  
                                       Mean   :100   Mean   :12762.7  
                                       3rd Qu.:100   3rd Qu.:18122.3  
                                       Max.   :100   Max.   :41395.4  

1.8 Computing Huff model

huff_share <- huff.shares(
  interaction_matrix, "VILLCODE", "Store_CD", 
  "Store_size", "distance")

1.8.1 Trade area analysis for a single store

First, let us select market share data of a store by usign Store_CD field (i.e. TG).

market_share <- huff_share %>%
  filter(Store_CD == "BT")

Next, trade areas of the selected will be derived by using the code chun below.

Trade_area <- village %>%
  left_join(market_share)

Joining with `by = join_by(VILLCODE)`

Before plotting the trade area map, code chunk below is used to extract the selected store from stores sf data.frame.

selected_store <- stores %>%
  filter(Store_CD == "BT")

Now, use the code chunk below to plot the trade areas of the selected store.

tm_shape(Trade_area) +
  tm_polygons(fill = "p_ij",
              fill.scale = tm_scale_intervals(
                style = "quantile", 
                n = 10,
                values = "brewer.blues"),
              fill.legend = tm_legend(
                title = "Patronage Probability")) +
  tm_shape(selected_store) +
    tm_dots(size = 0.3,
            fill = "red",
            col = "black") +
  tm_title("Distribution of Patronage Probability of Store BT by village") +
  tm_layout(frame = TRUE) +
  tm_borders(fill_alpha = 0.5) +
  tm_compass(type="8star", size = 2) +
  tm_scalebar() +
  tm_grid(alpha =0.2) 

1.8.2 Store sales analysis

One of the usage of Huff model is to estimate the store sales by using village population.

Firstly, let us append P_CNT field from village sf.data.frame to huff_share data.frame by using the code chunk below.

huff_share_pp <- huff_share %>%
  left_join(village %>% 
      st_drop_geometry() %>% 
      select(VILLCODE, P_CNT),
    by = "VILLCODE") %>%
  relocate(
    distance, .after = last_col())

Next, shares.total() of MCI package will be used to calculate the total sales by using the code chunk below.

huff_total <- shares.total(
  huff_share_pp, "VILLCODE", "Store_CD",
  "p_ij", "P_CNT")

Code chunk below is used to append the estimate store sales (Ej) onto stores sf data.frame.

store_sales_est <- stores %>%
  left_join(huff_total,
            by = join_by(
              "Store_CD" == "suppliers_single"))

Now, we can visualise how well the estimate sales as compare to the actual sales by using the code chunk below.

store_sales_nongeom <- store_sales_est %>% st_drop_geometry()

ggscatterstats(
  store_sales_nongeom,
  x = sum_E_j,
  y = Total_Sales)

`stat_xsidebin()` using `bins = 30`. Pick better value `binwidth`.
`stat_ysidebin()` using `bins = 30`. Pick better value `binwidth`.

huff_share_wider <- huff_share %>%
  select(VILLCODE, Store_CD, p_ij) %>%
    pivot_wider(names_from = Store_CD, 
                values_from = p_ij)